##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)

##########################################################################################

option_list <- list(
  make_option(c("--Gene"), type = "character") ,
  make_option(c("--stad_mut_rate_dgc"), type = "character") ,
  make_option(c("--stad_mut_rate_igc"), type = "character") ,
  make_option(c("--images_path"), type = "character")
)

if(1!=1){
  
  Gene <- "TP53"
  work_dir <- ""
  stad_mut_rate_dgc <- paste(work_dir,"/OncoSG_STAD_DGC.mut_rate.tsv",sep="")
  stad_mut_rate_igc <- paste(work_dir,"/OncoSG_STAD_IGC.mut_rate.tsv",sep="")
  images_path <- paste(work_dir,"/mutRatePlot",sep="")
  
}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

Gene <- opt$Gene
stad_mut_rate_dgc <- opt$stad_mut_rate_dgc
stad_mut_rate_igc <- opt$stad_mut_rate_igc
images_path <- opt$images_path


###########################################################################################
stad_mut_rate_dgc <- data.frame(fread(stad_mut_rate_dgc))
stad_mut_rate_dgc$SampleNum <- 33
stad_mut_rate_dgc$Class <- "DGC"
stad_mut_rate_igc <- data.frame(fread(stad_mut_rate_igc))
stad_mut_rate_igc$SampleNum <- 63
stad_mut_rate_igc$Class <- "IGC"
stad_mut_rate_dgc <- subset(stad_mut_rate_dgc , gene == Gene)
stad_mut_rate_igc <- subset(stad_mut_rate_igc , gene == Gene)
mut_rate_gene_file <- rbind(stad_mut_rate_dgc,stad_mut_rate_igc)

col <- c(
  rgb(234,106,79,alpha=255,maxColorValue=255),
  rgb(203,24,30,alpha=255,maxColorValue=255)
)

names(col) <- c("IGC" , "DGC")

###########################################################################################

dat_plot <- mut_rate_gene_file
dat_plot$Class <- factor( dat_plot$Class , levels = c("IGC" , "DGC") , order = T )
dat_plot$value_text <- paste0( round(dat_plot$rate , 2) * 100 , "%") 

###########################################################################################
## 计算P值
dat_plot$p.value = ""
dat_plot$p_text = ""

for(geneN in unique(dat_plot$gene)){
  
  print(geneN)
  
  tmp_1 <- subset( dat_plot , gene == geneN & Class %in% c("IGC") )
  tmp_2 <- subset( dat_plot , gene == geneN & Class %in% c("DGC") )
  
  if(nrow(tmp_1)==0){
    tmp_1 <- tmp_2
    tmp_1$MutNum <- 0
    tmp_1$MutRate <- 0
    tmp_1$value_text <- ""
  }
  
  if(nrow(tmp_2)==0){
    tmp_2 <- tmp_1
    tmp_2$MutNum <- 0
    tmp_2$MutRate <- 0
    tmp_2$value_text <- ""
  }
  
  tmp <- rbind( tmp_1 , tmp_2 )
  
  tmp_fisher <- matrix(c(tmp$freq , tmp$SampleNum - tmp$freq) , ncol = 2)
  p <- fisher.test(tmp_fisher)$p.value
  
  dat_plot[dat_plot$gene == geneN & dat_plot$Class %in% c("IGC" , "DGC"),"p.value"] <- p
  p_text <- paste0( "DGC vs IGC \np = " , format(as.numeric(p) , scientific = TRUE , digits = 3) )
  dat_plot[dat_plot$gene == geneN & dat_plot$Class %in% c("IGC" , "DGC"),"p_text"] <- p_text
}

dat_plot$gene <- factor( dat_plot$gene , 
                         levels = unique(dat_plot[order(dat_plot$freq , decreasing=T),"gene"]) , order = T)

###########################################################################################

dat_plot <- subset( dat_plot , gene == Gene )

dat_plot2 <- dat_plot
dat_plot2$Type <- ""
dat_plot2$Type[dat_plot2$Class %in% c("IGC" , "DGC")] <- "IGC vs DGC"
dat_plot2$p_text <- paste0( "p = " , format(as.numeric(dat_plot2$p.value) , scientific = TRUE , digits = 3))
dat_plot2$Type <- factor( dat_plot2$Type , levels = "IGC vs DGC" )

###########################################################################################
## IGC vs DGC
plot <- ggplot( data = dat_plot2[dat_plot2$Type=="IGC vs DGC",] , aes( x = Class , y = rate , fill = Class ))+
  geom_bar(position = "stack", stat = "identity") + 
  theme_bw()+
  labs(x="",y="Mutation Rate")+
  facet_wrap(~Type,scales="free")+
  theme(panel.grid = element_blank())+
  scale_fill_manual(values=col) +
  ylim(0,1.05)+
  geom_text(aes(label=p_text , y = 1 ,x = 1.5),size=3)+
  geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3 , fontface='bold' , color="black")+
  theme(panel.background = element_blank(),#设置背影为白色#清除网格线
        legend.position ='right',
        legend.title = element_blank() ,
        panel.grid.major=element_line(colour=NA),
        plot.title = element_text(size = 12,color="black",face='bold'),
        legend.text = element_text(size = 12,color="black",face='bold'),
        axis.text.y = element_text(size = 7,color="black",face='bold'),
        axis.title.x = element_text(size = 12,color="black",face='bold'),
        axis.title.y = element_text(size = 12,color="black",face='bold'),
        strip.text.x = element_text(size = 7 , face = 'bold'),
        axis.ticks.x = element_blank(),
        axis.text.x = element_text(size = 8,color="black",face='bold') ,
        )

out_name <- paste0( images_path , "/OncoSG_MutRate_" , Gene , ".IGC_DGC.pdf" )
ggsave( out_name , plot , width = 3.5 , height = 5 )
